/**
 * @file QuadraticInterpolation.h
 * @brief The implementation of the Quadratic interpolation operator.
 * @author XDDDDD
 * @version 
 * @date 2021-06-10
 */

#ifndef __PAISLEYPARK__QUADRATICINERPOLATION_H__
#define __PAISLEYPARK__QUADRATICINERPOLATION_H__

#include "InterpolationOperator.h"

template <size_t DIM>
class QuadraticInterpolation : public InterpolationOperator<DIM> {
public:
	void operator()(MultigridVector<DIM>& _v, size_t _M) override;
};

template <size_t DIM>
void QuadraticInterpolation<DIM>::operator()(MultigridVector<DIM>& _v, size_t _M) {
		if(_M == 0 || _M >= _v.get_N()) {
		std::cout << "Invalid _M" << std::endl;
		std::exit(-1);
	}
	std::vector<size_t> coord_2h, coord_h;
	coord_2h.clear();
	coord_2h.resize(DIM);
	coord_h.clear();
	coord_h.resize(DIM);
	std::vector<size_t> coord;
	coord.clear();
	coord.resize(DIM);
	int d[DIM] = {0};
	size_t pole = 0;
	size_t num_1d = (1 << _M) + 1;
	double coef = 1.0;
	double w = 0;
	bool is_jump;
	bool is_mid;
	bool is_out;
	for(size_t i = 0; i != DIM; i++) {
		coord_2h[i] = 1;
	}
	//clear.
	_v.cleanpart(_M + 1);
	//Go through all the nodes in the coareser grid.
	while(pole != DIM) {
		while(pole != DIM) {
			if(_M == 1) {
				pole = 0;
				is_mid = true;
				coef = 1.0;
				for(size_t i = 0; i != DIM; i++) {
					if(d[i] != 0) {
						is_mid = false;
						coef = coef*0.5;
					}
					coord_h[i] = 2 + d[i];
					coord_2h[i] = 1;
				}
				if(is_mid)
					_v(coord_h, _M + 1) = _v(coord_2h, _M);
				else {
					_v(coord_h, _M + 1) += coef*6.0*0.25*_v(coord_2h, _M);
					for(size_t i = 0; i != DIM; i++) {
						if(d[i] != 0)
							coord_2h[i] = 0;
						_v(coord_h, _M + 1) -= coef*0.25*_v(coord_2h, _M);
					}
					for(size_t i = 0; i != DIM; i++) {
						if(d[i] != 0)
							coord_2h[i] = num_1d*2 - 2;
						_v(coord_h, _M + 1) -= coef*0.25*_v(coord_2h, _M);
					}
				}
			}
			else {
				pole = 0;
				is_mid = true;
				is_out = false;
				coef = 1.0;
				for(size_t i = 0; i != DIM; i++) {
					if(d[i] != 0) {
						is_mid = false;
						coef = coef*0.5;
					}
					coord_h[i] = 2*coord_2h[i] + 3*d[i];
				}
				if(is_mid)
					_v(coord_h, _M + 1) = _v(coord_2h, _M);
				else {
					for(size_t i = 0; i != DIM; i++) {
						if(!(coord_h[i] > 0 && coord_h[i] < 2*num_1d - 2)) {
							is_out = true;
						}
					}
					if(is_out) {
						for(size_t i = 0; i != DIM; i++) {
							coord_h[i] = coord_2h[i]*2 + d[i];
							coord_2h[i] = coord_2h[i] - d[i];
						}
						_v(coord_h, _M + 1) -= coef*0.125*_v(coord_2h, _M);
						for(size_t i = 0; i != DIM; i++) {
							coord_2h[i] = coord_2h[i] + 2*d[i];
						}
						_v(coord_h, _M + 1) += coef*0.25*3*_v(coord_2h, _M);
						for(size_t i = 0; i != DIM; i++) {
							coord_2h[i] = coord_2h[i] - d[i];
						}
						_v(coord_h, _M + 1) += coef*0.25*6*_v(coord_2h, _M);
					}
					else {
						for(size_t i = 0; i != DIM; i++) {
							coord_h[i] = coord_2h[i]*2 + 3*d[i];
						}
						_v(coord_h, _M + 1) -= coef*0.125*_v(coord_2h, _M);
						for(size_t i = 0; i != DIM; i++) {
							coord_h[i] = coord_2h[i]*2 + d[i];
						}
						_v(coord_h, _M + 1) += coef*9.0*0.125*_v(coord_2h, _M);
					}
				}
			}
			d[pole] = -2*d[pole] - 1;
			while(d[pole] == -3) {
				d[pole] = 0;
				pole++;
				if(pole == DIM)
					break;
				d[pole] = -2*d[pole] - 1;
			}
		}
		pole = 0;
		if(_M == 1)
			pole++;
		else {
			coord_2h[pole]++;
			while(coord_2h[pole] == num_1d - 1) {
				coord_2h[pole] = 1;
				pole++;
				if(pole == DIM)
					break;
				coord_2h[pole]++;
			}
		}
	}
};

#else
//Do nothing.
#endif
